#define _GNU_SOURCE
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#include "netpbm.h"

/* Reemplazo de las constantes por define de preprocesamiento:
    Nos da una mejora sustancia con respecto a los stalled-cycles
    y mejora ligeramente el tiempo de ejecución */

#define N 500
#define SOURCE_TEMP 5000.0f
#define BOUNDARY_TEMP 2500.0f
#define MIN_DELTA 0.05f
#define MAX_ITERATIONS 10000

static unsigned int idx(unsigned int x, unsigned int y, unsigned int stride) {
    return y * stride + x;
}


static void init(unsigned int source_x, unsigned int source_y, float * matrix) {
    // init
    memset(matrix, 0, N * N * sizeof(float));

    // place source
    matrix[idx(source_x, source_y, N)] = SOURCE_TEMP;

    // fill borders
    for (unsigned int x = 0; x < N; ++x) {
        matrix[idx(x, 0,   N)] = BOUNDARY_TEMP;
        matrix[idx(x, N-1, N)] = BOUNDARY_TEMP;
    }
    for (unsigned int y = 0; y < N; ++y) {
        matrix[idx(0,   y, N)] = BOUNDARY_TEMP;
        matrix[idx(N-1, y, N)] = BOUNDARY_TEMP;
    }
}


static void step(unsigned int source_x, unsigned int source_y, const float * current, float * next) {
    for (unsigned int y = 1; y < N-1; ++y) {
        for (unsigned int x = 1; x < N-1; ++x) {
            if ((y == source_y) && (x == source_x)) {
                continue;
            }
            next[idx(x, y, N)] = (current[idx(x, y-1, N)] +
                                  current[idx(x-1, y, N)] +
                                  current[idx(x+1, y, N)] +
                                  current[idx(x, y+1, N)]) / 4.0f;
        }
    }
}


static float diff(const float * current, const float * next) {
    float maxdiff = 0.0f;
    for (unsigned int y = 1; y < N-1; ++y) {
        for (unsigned int x = 1; x < N-5; x+=5) {
            maxdiff = fmaxf(maxdiff, fabsf(next[idx(x, y, N)] - current[idx(x, y, N)]));
            maxdiff = fmaxf(maxdiff, fabsf(next[idx(x+1, y, N)] - current[idx(x+1, y, N)]));
            maxdiff = fmaxf(maxdiff, fabsf(next[idx(x+2, y, N)] - current[idx(x+2, y, N)]));
            maxdiff = fmaxf(maxdiff, fabsf(next[idx(x+3, y, N)] - current[idx(x+3, y, N)]));
            maxdiff = fmaxf(maxdiff, fabsf(next[idx(x+4, y, N)] - current[idx(x+4, y, N)]));
        }
        maxdiff = fmaxf(maxdiff, fabsf(next[idx(N-5, y, N)] - current[idx(N-5, y, N)]));
        maxdiff = fmaxf(maxdiff, fabsf(next[idx(N-4, y, N)] - current[idx(N-4, y, N)]));
        maxdiff = fmaxf(maxdiff, fabsf(next[idx(N-3, y, N)] - current[idx(N-3, y, N)]));
        maxdiff = fmaxf(maxdiff, fabsf(next[idx(N-2, y, N)] - current[idx(N-2, y, N)]));
    }
    return maxdiff;
}


int main() {
    size_t array_size = N * N * sizeof(float);

    float * current = malloc(array_size);
    float * next = malloc(array_size);

    srand(time(NULL));
    unsigned int source_x = rand() % (N-2) + 1;
    unsigned int source_y = rand() % (N-2) + 1;
    printf("Heat source at (%u, %u)\n", source_x, source_y);

    init(source_x, source_y, current);
    memcpy(next, current, array_size);

    float t_diff = SOURCE_TEMP;
    for (unsigned int it = 0; (it < MAX_ITERATIONS) && (t_diff > MIN_DELTA); ++it) {
        step(source_x, source_y, current, next);
        t_diff = diff(current, next);
        printf("%u: %f\n", it, t_diff);

        float * swap = current;
        current = next;
        next = swap;
    }

    char * filename;
    asprintf(&filename, "heat.ppm");
    writePPMbinaryImage(filename, current, N, N, N, 0.0f, fmaxf(SOURCE_TEMP, BOUNDARY_TEMP));
    free(filename);

    free(current);
    free(next);

    return 0;
}
